clear all
close all

% =====================================================================
% ====  "Risk-Sharing in Village Economies Revisited", by Bold and Broer

% ====   This programme solves the limited commitment village economy for
% both the standard specification (individual deviations) and the
% alternative specification (coalitional deviations)

% This version to be used for all vilages and cases apart from Shirapur with colitional
% deviations

% ====    Written by Tessa Bold and Tobias Broer

% ====    This version Septermber 2021
% =====================================================================


filetosave='some_name';            % name of workspace saved
Tessa1Tobi2=1;                     % sets paths depending on user below
if Tessa1Tobi2==1
currentfolder='C:\Users\tbold\Dropbox\Tessa and Tobi\JEEA Replication';    
elseif Tessa1Tobi2==2
currentfolder='C:\Users\tbroer\Dropbox\Research\Current_Projects\tessa\JEEA Replication';
end

%====================================================
% select parameters
% ===================================================

village=2;                         % village indicator: 1 - Aurepalle; 2-Kanzara; 3-Shirapur
coaldev=1;                         % set to 1 for coalitional deviations, to 2 for self-insurance (SI), to 0 for individual deviations (ID), to 3 for full-insurance with coalitional deviations
unconditional=0;                   % benchmark=0; set this to 1 for unconditional income data, to 0 for income data conditional on time fixed effets / demeaned data period-by-period
fixedeffect=1;                     % benchmark=1; set this to 1 for an income process with fixed effects
one_per_aut=0;                     % benchmark=0; set this to 1 if you want deviations to imply one period of autarky
ytaxrate=0;                        % benchmark=0; linear tax on income
scatter = 0;                       % benchmark=0; set to 1 to just draw a scatter plot
Rov_sb = 1;                        % benchmark=0; set this to 1 if you want the outside option for the Rov in the case of coalitional deviations to equal the average value from the solution to the next smaller village (second-best), rather than
% just the utility from consuming average ROV income (first-best)
estrho1=0;                         % benchmark=0; set this to 1 if you would like to estimate persistence, default is 0
rhogrid=0;                         % benchmark=0; set to 1 if you want to have also a grid for CRRA rho - only to show non-identification in Aurepalle (village=1)
Nphi=1;                            % number of borrowing limits in SI economy - not used currently
server=0;                          % set to 1 to run on server; o wise path for workstations is used
buildup=0;                         % benchmark=0; set to 1  to calculate ID model for subset of group sizes, not just whole village
Nincproc=1;                       % benchmark=10; number of income processes to look at for income measurement error
if estrho1==0
    maxincerr=2/3;                     % maximum income measurment error in multiples of income variance (in levels)
elseif estrho1==1
    maxincerr=0;
end
model_var_dy=1;                    % set to 1 if you want to use the model-implied variance of dlog(y) to specify maximum meas error, otherwise var(dlog(y)) from data is used, slightly different as we target var(log(y)), not var(dlog(y())
T=6200; Tinit=201;                 % number of simulation periods in total, T-Tinit+1 equals number of initial periods to discard
maxcerr=0.9;                       % max cons meas error in log levels
if scatter==0
    N_cerr=30;                         % number of support points for cons meas error
elseif scatter==1
    N_cerr=0;
end
if coaldev==0
    convlength=0;                  % the conv criterion must be fulfilled at least in convlength consecutive periods to avoid jumping 'in'
else
    convlength=15;
end
if village==1
    vss_high_inddev=34-1;          % the number of total village members is vss_high_inddev+1 under individual deviations
    vss_high_coaldev=34-1;
elseif village==2;
    vss_high_inddev=37-1;
    vss_high_coaldev=37-1;
elseif village==3
    vss_high_inddev=31-1;
    vss_high_coaldev=31-1;
end
vssbuildupvec_inddev=[1:7,10,15,vss_high_inddev]; % for buildup==1, this is the vector of village sizes for which the solution is calculated in the individual deviations case
if server==1
   
else
    if Tessa1Tobi2==2
        outputpath=[currentfolder '\Programmes/Main Results and Robustness']
        programpath=[currentfolder '\Programmes']
    elseif Tessa1Tobi2==1
        outputpath=[currentfolder '\Programmes/Main Results and Robustness']
        programpath=[currentfolder '\Programmes']
    end
end
cd(programpath)
if buildup==0
    if rhogrid==0 && scatter == 0
        betavec0=[0.5:0.1:0.7,0.8:0.01:0.98];
    elseif scatter ==1
        betavec0=[0.83];
    elseif rhogrid==1
        betavec0=[0.83,0.85,0.87,0.89,0.91,0.93,0.95,0.97];
    end
elseif buildup==1
    betavec0=[0.86,0.91]; 
end
if rhogrid==0
    if scatter==1
        betavec1=[0.97];
    else
        betavec1=[0.84:0.02:0.88,0.89:0.01:0.99];
    end
elseif rhogrid==1
    betavec1=[0.84,0.86,0.88,0.90,0.92,0.94,0.96,0.98];
end
betavec3=[0.96:0.0025:0.99];
if coaldev<2 || coaldev==3
    rhovec=[1];
elseif coaldev==2
    rhovec=ones(1,Nphi);
end

if rhogrid==0
    gapcritset=0.0001;              %convergence criterion for change in values
elseif rhogrid==1
    gapcritset=0.001;
end

rov_average=1;                      % set to 1 to specify the rest of village in terms of average
% consumption
momentclevinlog=1;                  % set this to 1 if you want to calculate the moments for levels of
% consumption based on log consumption
maxitgapset=400;                    % maximum number of iterations
Pun_fac=0.0;                        % consumption penalty in the first period after deviation equal
% to caut*Pun_fac
% grid of time-varying planner weights
lambda=[0.05:0.001:0.95]';
m=length(lambda);
if floor((m+1)/2)~=((m+1)/2)
    disp('care: lambda has to have uneven no of points in first dimension')
    stop
end
S_KEEP=nan(1000*3,2,vss_high_inddev);               % pre-allocate the matrices for coalitional deviations, so that you don't run simulations vss times
S_rov_KEEP=nan(100,1,vss_high_inddev);              % this one is for a maximum of vss=19 for the CD
P_rov_KEEP=nan(100,1000,vss_high_inddev);
PHI_KEEP=nan(m,1000*3,2,vss_high_inddev);
IDIO_DIST_ROV_KEEP=nan(1000*3,3,vss_high_inddev);

% ====================================================
% Specify Income process for individual and Rest of Village (rov)
% ====================================================
% data moments for log income from ICRISAT data - all villages
% first col: Aurepalle, second: Kanzara, third: Shirapur
if fixedeffect==1
    if unconditional==0
        datafile='armoments_fe0';
    elseif unconditional==1
        datafile='armoments_fe1';
    end
elseif fixedeffect==0
    datafile='armoments';
end
if server==0
    cd([currentfolder '\Data\Output']) % need to change this according to user
end
eval(['load ' datafile '.out;']);
eval(['covyylag=' datafile '(1,:);']);
eval(['vardy=' datafile '(2,:);']);
eval(['vary=' datafile '(3,:);']);
if model_var_dy==1
    vardy=[0.231,0.125, 0.269];
end
cd(programpath)
if maxincerr>0
    if Nincproc>1
        varnu_max=maxincerr*vardy/2;                % nu is measurement error in log levels
    else
        varnu_max=zeros(1,3);
    end
    rho_min=covyylag(village)/vary(village);                      % AR(1) coefficient w/o meas error
    rho_max=covyylag(village)/(vary(village)-varnu_max(village)); % AR(1) coefficient w max meas error
    rho1vec=linspace(rho_min,rho_max,Nincproc);
    varnu=max(vary(village)-covyylag(village)./rho1vec,0); % max prevents the first entry to be -e17
    vareps=(vary(village)-varnu).*(1-rho1vec.^2);
else
    rho1_data=covyylag(village)/vary(village);
    vareps_data=vary(village).*(1-rho1_data.^2);
    rho1vec=[-0.6,-0.4,-0.2,0,0.2,0.4,0.6];
    vareps=vary(village).*(1-rho1vec.^2);
    varnu=zeros(1,length(vareps));
    Nincproc=length(vareps);
end

% ===================================
% Iterate over income processes
% ===================================
disp('now solving for village,coaldev, unconditional, Nincproc, fixedeffect one_per_aut Rov_sb - press any key to continue')
disp([village coaldev unconditional Nincproc fixedeffect one_per_aut Rov_sb])
pause
indicator1=zeros(1,size(rhovec,2));
for incproc=1:Nincproc
    if rhogrid==0 && buildup==0 && maxincerr>0 && Rov_sb==0 && one_per_aut==1 && scatter == 0 %
        filetosaveest=['Momentsest' int2str(village) int2str(coaldev) int2str(incproc) int2str(unconditional) int2str(fixedeffect) int2str(one_per_aut)];
    elseif maxincerr==0
        filetosaveest=['EstrhoMomentsest' int2str(village) int2str(coaldev) int2str(incproc) int2str(unconditional) int2str(fixedeffect) int2str(one_per_aut)];
    elseif Rov_sb==1
        filetosaveest=['SBMomentsest' int2str(village) int2str(coaldev) int2str(incproc) int2str(unconditional) int2str(fixedeffect) int2str(one_per_aut)];
    elseif rhogrid==1
        filetosaveest=['GridMomentsest' int2str(village) int2str(coaldev) int2str(incproc) int2str(unconditional) int2str(fixedeffect) int2str(one_per_aut) ];
    elseif buildup==1
        filetosaveest=['BuildupMomentsest' int2str(village) int2str(coaldev) int2str(incproc) int2str(unconditional) int2str(fixedeffect) int2str(one_per_aut)];
    elseif scatter==1
        filetosaveest=['ScatterMomentsest' int2str(village) int2str(coaldev) int2str(incproc) int2str(unconditional) int2str(fixedeffect) int2str(one_per_aut)];
    elseif one_per_aut==0   % this is the benchmark option
        filetosaveest=['No_1_per_aut_Momentsest' int2str(village) int2str(coaldev) int2str(incproc) int2str(unconditional) int2str(fixedeffect) int2str(one_per_aut)];
    end
    
    % specify a markov process for individual income
    [incomevals, Q1] = rouwenhorst(rho1vec(incproc),vareps(incproc)^0.5,3);
    incomevals=exp(incomevals');
    QQQ=Q1^1000;
    SS=QQQ(1,:)';
    incomevals=incomevals/(SS'*incomevals);
    cumQ=cumsum(Q1,2);
    incomevals=(1-ytaxrate)*incomevals+ytaxrate; % tax as in Krueger and Perri JET 2011
    
    % set minimum and maximum village size
    if coaldev==1
        if rhogrid==0
            vss_set=[1:7];% prior solutions have shown that 1-7 is the relevant set. Otherwise use, for example, [1:12,17,22,27,vss_high_coaldev]
        elseif rhogrid==1
            vss_set=[1:5]
        end
        vss_high=vss_set(end);
    elseif coaldev==0
        if buildup==0
            vss_set=[vss_high_inddev];
            lagind=1;
            vss_high=vss_set(end);
        elseif buildup==1
            vss_set=vssbuildupvec_inddev;
            lagind=1;
            vss_high=vss_set(end);
        end
    elseif coaldev==2;
        vss_set=vss_high_inddev;
        vss_high=vss_set(end);
    elseif coaldev==3
        vss_set=[1:10];
        vss_high=10;
    end
    
    % ===================================
    % iterate over discount factors
    % ===================================
    if coaldev==0
        betavec=betavec0;
    elseif coaldev==1
        betavec=betavec1;
    elseif coaldev==2
        r=0.04;
        if scatter==0
            betamax=1/(1+r)-0.015;
            betavec=linspace(0.8,betamax,20);
        elseif scatter ==1
            betavec=0.92;
        end
        Z=incomevals;
        phinat=min(Z)/(r-1);
        if Nphi>1
            phivec=linspace(phinat,0,Nphi);
        else
            phivec=0;
        end
    elseif  coaldev==3
        betavec=betavec3;
    end
    Momentsest=nan(35,N_cerr+1,size(betavec,2),1,vss_high_inddev);
    for Qbeta=1:size(betavec,2)
        beta=betavec(Qbeta);
        if rhogrid==1 % this is only for the identification graph
            if coaldev==0
                if beta==0.83
                    rhovec=1.468461538461538
                elseif beta==0.85
                    rhovec=1.316666666666667
                elseif beta==0.87
                    rhovec=1.16
                elseif beta==0.89
                    rhovec=1
                elseif beta==0.91
                    rhovec=0.8340
                elseif beta==0.93
                    rhovec=0.6611
                elseif beta==0.95
                    rhovec=0.4808
                elseif beta==0.97
                    rhovec=0.2912
                end
            elseif coaldev==1
                if beta == 0.98
                    rhovec=0.366222222222222
                elseif beta==0.84
                    rhovec=2.389
                elseif beta==0.86
                    rhovec=2.11
                elseif beta==0.88
                    rhovec=1.821
                elseif beta==0.90
                    rhovec=1.5515
                elseif beta==0.92
                    rhovec=1.28718
                elseif beta==0.94
                    rhovec=1
                elseif beta==0.96
                    rhovec=0.6935
                else
                end
            end
        end
        
        % ===================================
        % iterate over CRRA
        % ===================================
        for Qrho=1:size(rhovec,2)
            rho=rhovec(Qrho)  ;
            if coaldev==1
                lagind=[];
            else
                lagind=1;
            end
            P_LAG=nan(1000,1000,length(vss_set));Y_lag=nan(1000,1000,length(vss_set));
            waut_lag=[]; waut_avg_lag=[]; waut_rov_sb=[]; waut_rov=[];waut_temp=[];
            WAUT_LAG=nan(1000,1000,length(vss_set));
            STABLE=[];
            if coaldev==3
                Stable=nan(length(vss_set));
                stable=[]; countstable=0;
                clear Delta
            end
            for vsscount=1:length(vss_set)       % iterate over village size
                
                vss=vss_set(vsscount);
                vss_max=max(vss_set);
                if vsscount>1 & coaldev==1
                    lagind=lagind+(vss-(vss_set(vsscount-1)+1));
                end
                disp('now solving for parameters')
                disp([Qbeta,Qrho])
                disp('vss')
                disp(vss)
                % Aiyagari economy
                if coaldev==2
                    sav_problem
                else
                    % this is the scaling factor: if you want RoV consumption as average consumption per HH, scale rov variables
                    % by vss
                    if rov_average==1
                        ScaleRov=vss;
                    else
                        ScaleRov=1;
                    end
                    csim=[];  cc =[]; eeewaut =[]; eewaut =[]; ewaut =[]; waut =[]; waut_check=[]; waut_check_ind=[]; rr =[]; cshare=[]; raw =[]; cshare =[]; dcshare =[]; dcshare_raw =[]; S_rov =[]; P_rov =[];
                    S_rov_new =[]; P_rovnew =[]; S =[]; P =[]; caut =[]; c =[]; waut =[]; gammagrid =[]; gammar =[]; phi =[]; phisol =[]; csol =[];csolold=[]; cfirstbest =[];
                    wfirstbest =[]; ewfirstbest =[]; wsbnewsol=[];
                    Y =[]; vs=[];  idio_dist_autnew=[];idio_dist_aut_new=[];normperf=[];w1=[];indSj=[]; indSjj=[]; Ssim=[]; aggvillincome=[];aggincgrowth=[];aggvillincgrowth=[];
                    S_rov_aut=[]; S_rov_sb=[]; S_rov_aut_all=[];S_rov_aut_all2=[];
                    w =[];  phinew =[]; Yagg =[]; wint =[]; wsb =[]; ewsb =[]; csol =[]; phisol =[]; nonconv =[];   gammarfun =[];  x0 =[]; index =[]; ewsb =[]; bind=[];
                    waut=[];
                    
                    % ===================================
                    % build the state space of ROV income
                    % ===================================
                    idio_dist=[];
                    idio_dist_aut=[];
                    idio_dist_aut_all=[];idio_dist_aut_all2=[];
                    idio_dist_sb=[];
                    
                    for k1=0:vss
                        for k2=0:vss
                            for k3=0:vss
                                II=zeros(1,3);
                                if k1+k2+k3==vss
                                    % the vector of aggregate incomes from
                                    % permutations of vss households
                                    S_rov=[S_rov; k1*incomevals(1,1) + k2*incomevals(2,1) + k3*incomevals(3,1)];
                                    % the corresponding distribution of
                                    % individual income values
                                    idio_dist=[idio_dist; k1 k2 k3];
                                    % this computes the 'best' distribution of incomes
                                    % of size vss-lagind
                                    if vss>1
                                        II=zeros(1,3);
                                        temp=[k1 k2 k3];
                                        % this takes away the lowest incomes to
                                        % derive the distribution of incomes of the
                                        % 'best' coalition lagind individuals
                                        % This coalition is one smaller than the
                                        % largest sustainable coalition because it
                                        % will consist of 1 agent and m-1 guys from
                                        % the largest sustainable coalition of size
                                        % m
                                        for jj=1:lagind
                                            II(min(find(temp>0)))=II(min(find(temp>0)))+1;
                                            temp=[k1 k2 k3]-II;
                                        end
                                        idio_dist_aut=[idio_dist_aut; [k1 k2 k3]-II];
                                        S_rov_aut=[S_rov_aut; ([k1 k2 k3]-II)*incomevals];
                                        
                                        II=zeros(1,3);
                                        temp=[k1 k2 k3];
                                        % this takes away the lowest incomes to
                                        % derive the distribution of incomes of the
                                        % 'best' coalition
                                        % consisting of m individuals from only the
                                        % rest of village
                                        for jj=1:lagind-1
                                            II(min(find(temp>0)))=II(min(find(temp>0)))+1;
                                            temp=[k1 k2 k3]-II;
                                        end
                                        idio_dist_sb=[idio_dist_sb; [k1 k2 k3]-II];
                                        S_rov_sb=[S_rov_sb; ([k1 k2 k3]-II)*incomevals];
                                    end
                                end
                            end
                        end
                    end
                    
                    % resort S_rov in ascending order
                    [S_rov,S_rovind]=sort(S_rov);
                    idio_dist=idio_dist(S_rovind,:);
                    
                    if vss>1
                        idio_dist_aut=idio_dist_aut(S_rovind,:);
                        S_rov_aut=S_rov_aut(S_rovind);
                        idio_dist_sb=idio_dist_sb(S_rovind,:);
                        S_rov_sb=S_rov_sb(S_rovind);
                    end
                    
                    %%%%and create S_rov_aut vectors that drops lagind individuals from
                    %%%%the rov vector.
                    if vss>1
                        for k=1:length(STABLE)
                            if STABLE(k)==1
                                for g=1:length(idio_dist)
                                    index=0;
                                    for j1=0:vss-k
                                        for j2=0:vss-k
                                            for j3=0:vss-k
                                                if j1+j2+j3==vss-k && (idio_dist(g,1)>=j1 & idio_dist(g,2)>=j2 & idio_dist(g,3)>=j3)
                                                    index=index+1;
                                                    temp=[j1 j2 j3];
                                                    idio_dist_aut_all(g,:,index,k)=[idio_dist(g,:)-temp];
                                                    S_rov_aut_all(g,:,index,k)=[(idio_dist(g,:)-temp)*incomevals];
                                                end
                                            end
                                        end
                                    end
                                end
                            end
                            
                        end
                    end
                    
                    % ===================================
                    % simulate aggregate income transitions
                    % ===================================
                    tic
                    T_trans=50000;
                    aggincind=zeros(length(S_rov),length(S_rov));
                    for jjj=1:length(S_rov)
                        jjj;
                        inccount=zeros(T_trans,vss,3);
                        agginc=zeros(T_trans,3);
                        for iii=find(idio_dist(jjj,:)>0)
                            for kkk=1:idio_dist(jjj,iii)
                                rng(jjj*100+iii*10+kkk,'twister');
                                shock=rand(T_trans,1);
                                rr=(find(shock<=cumQ(iii,1)));
                                inccount(rr,kkk,iii)=1;
                                rr=(find(cumQ(iii,1)<shock & cumQ(iii,2)>=shock));
                                inccount(rr,kkk,iii)=2;
                                rr=(find(cumQ(iii,2)<shock & cumQ(iii,3)>=shock));
                                inccount(rr,kkk,iii)=3;
                            end
                            agginc(:,iii)=sum(incomevals(inccount(:,1:idio_dist(jjj,iii),iii)),2);
                        end
                        agginc=sum(agginc,2);
                        for jjj2=1:length(S_rov)
                            qqq=(agginc<=S_rov(jjj2)+0.0000000001).*(agginc>=S_rov(jjj2)-0.0000000001);
                            aggincind(jjj,jjj2)=sum(qqq);
                        end
                        P_rov(jjj,:)=aggincind(jjj,:)/sum(aggincind(jjj,:));
                    end
                    disp('Simulation for agg income process takes')
                    toc
                    
                    %%scale everything correctly
                    if vss>1 && coaldev<2
                        S_rov_aut=S_rov_aut/(max(ScaleRov-lagind,1));
                        S_rov_sb=S_rov_sb/(max(ScaleRov-lagind+1,1));
                        for g=vss_set(find(vss-1>=vss_set))
                            if STABLE(g)==1
                                S_rov_aut_all(:,:,:,g)=S_rov_aut_all(:,:,:,g)/(g);
                            end
                        end
                    end
                    S_rov=S_rov/ScaleRov;
                    
                    % ===================================
                    % state space: combinations of indiv and village income
                    % ===================================
                    S=[kron(incomevals,ones(size(S_rov))),kron(ones(size(incomevals)),S_rov)];
                    idio_dist_rov=kron(ones(size(incomevals)),idio_dist);
                    idio_dist_aut=kron(ones(size(incomevals)),idio_dist_aut);
                    S_rov_aut=[kron(ones(size(incomevals)),S_rov_aut)];
                    S_rov_sb=[kron(ones(size(incomevals)),S_rov_sb)];
                    
                    if vss>1 && coaldev<2
                        for g=vss_set(find(vss-1>=vss_set))
                            if STABLE(g)==1
                                for k=1:size(S_rov_aut_all,3)
                                    S_rov_aut_all2(:,:,k,g)=[kron(ones(size(incomevals)),S_rov_aut_all(:,:,k,g))];
                                    idio_dist_aut_all2(:,:,k,g)=[kron(ones(size(incomevals)),idio_dist_aut_all(:,:,k,g))];
                                end
                            end
                        end
                        S_rov_aut_all=S_rov_aut_all2;
                        idio_dist_aut_all=idio_dist_aut_all2;
                        clear S_rov_aut_all2;
                        clear idio_dist_aut_all2;
                    end
                    
                    P=kron(Q1,P_rov);
                    
                    Y=S(:,1)+ScaleRov*S(:,2);			% possible aggregate income outcome
                    state=length(S);                    % number of states
                    agent=2;                            % number of 'agents' in the HH vs. RoV approximation equals 2
                    STATE(vss)=state;
                    
                    % ===================================
                    % Compute set of autarky values
                    % ===================================
                    % autarky values: coalitional deviations imply that the individual's
                    % outside option is given by utility of caut for 1 period, and then the
                    % expected value being in the 'best' insurance coalition of
                    % size vss-1, with most 'equal' weights, i.e.
                    % gamma=1/(vss-1) or, if this does not yield the autarky value corresponding to a coalition of vss-2, the value that is closest to it
                    for k=1:agent
                        caut(:,k)=S(:,k);
                    end
                    waut=zeros(state,agent);
                    waut_avg_lag=zeros(state,agent);
                    waut_sb=zeros(state,agent);
                    % this is just the PDV of utility from consumption of
                    % income caut
                    eeewaut(:,1)= inv(diag(ones(1,state),0)-beta.*P)*utility(rho,caut(:,1));
                    eeewaut(:,2)= inv(diag(ones(1,state),0)-beta.*P)*utility(rho,caut(:,2));
                    if coaldev<2
                        
                        % calculate the autarky value for ROV from average utility
                        % of its members in village of size vss
                        if vss>1 && Rov_sb==1
                            for j=1:state/3
                                % loop through individuals in the ROV
                                if lagind==1
                                    for jj=1:length(idio_dist(j,:))
                                        if idio_dist(j,jj)>0
                                            % distribution of incomes in ROV without one
                                            % individual
                                            idio_dist_lag=idio_dist(j,:);
                                            idio_dist_lag(1,jj)=idio_dist_lag(1,jj)-1;
                                            % make a temporary new 'state' of individual
                                            % and rov
                                            S_lag_temp(jj,2)=idio_dist_lag*incomevals;
                                            S_lag_temp(jj,1)=incomevals(jj);
                                            % find the corresponding state in a village of
                                            % size vss-1
                                            S_lag_2_temp=( ones(length(S_lag)/3,1)*idio_dist_lag==IDIO_DIST_ROV_KEEP(1:length(S_lag)/3,:,vss-1)) ;
                                            S_lag_2=sum(S_lag_2_temp,2)==3;
                                            S_lag_2=[S_lag_2;S_lag_2;S_lag_2];
                                            S_lag_1=(S_lag_temp(jj,1)==S_lag(:,1));
                                            S_lag_ind=find(S_lag_1.*S_lag_2==1);
                                            % construct utility of individual jj
                                            if one_per_aut==1
                                                waut_temp(jj,1)=utility(rho,S_lag_temp(jj,1)*(1-Pun_fac))+beta*P_lag(S_lag_ind,:)*waut_lag(:,1);
                                            else
                                                waut_temp(jj,1)=waut_lag(S_lag_ind,1);
                                            end
                                        else
                                            waut_temp(jj,1)=0;
                                        end
                                    end
                                    
                                    % utility of individuals is average of its members
                                    waut_rov(j,1)=idio_dist(j,:)*waut_temp/sum(idio_dist(j,:));
                                    
                                    % for the second best outside options, we don't look at cases with 'holes' in the sustainable
                                    % coalitions
                                elseif lagind>1
                                    waut_rov(j,1)=100000;
                                end
                            end
                            waut_rov_sb=[waut_rov;waut_rov;waut_rov];
                        end
                        
                        % in first iteration, or for individual deviations, waut is just consumption of income
                        if vss==1 || coaldev==0
                            for j=1:state
                                if one_per_aut==1
                                    waut(j,1)=utility(rho,caut(j,1)*(1-Pun_fac))+beta*P(j,:)*eeewaut(:,1);
                                    % autarky values for the rest of village are always the same - perfect
                                    % risk-sharing among the rest of village
                                    if Rov_sb==0 || vss==1
                                        waut(j,2)=utility(rho,caut(j,2)*(1-Pun_fac))+beta*P(j,:)*eeewaut(:,2);
                                    else
                                        waut(j,2)=utility(rho,caut(j,2)*(1-Pun_fac))+beta*P(j,:)*waut_rov_sb;
                                    end
                                else
                                    waut(j,1)=eeewaut(j,1);
                                    % autarky values for the rest of village are always the same - perfect
                                    % risk-sharing among the rest of village
                                    if Rov_sb==0 || vss==1
                                        waut(j,2)=eeewaut(j,2);
                                    else
                                        waut(j,2)=waut_rov_sb(j,1);
                                    end
                                    
                                end
                            end
                        else
                            % this calculates autarky values for the coalitional
                            % deviations case
                            %%%if you were going to use the second-best constrained
                            %%value, you first need to create the average lagged
                            %%utility weighted according to income of agent 1 and rest
                            %%of the village from the previous iteration
                            %1) look for the entry in S_rov_sb, which has the same
                            %average/aggregate income (1*S_lag(:,1)+(vss-lagind)*S_lag(:,2))
                            %2) this is usually more than one entry, since S_lag
                            %differs by which income agent 1 has and which income
                            %is rest of village. Here, we don't care about that
                            %anymore, only about the aggregate. Therefore we also
                            %need to average utility (hardly matters for which
                            %point you do this, to be really precise should take
                            %half form both), but important is to average utility
                            %as (waut_lag(ind_sb,1)+(vss-lagind)*waut_lag(ind_sb,2))/(vss-lagind+1)
                           
                            
                            for j=1:state
                                % indicator for the 'autarky' state in the state
                                % space of the previous iteration S_lag
                                S_lag_2=(abs(S_rov_aut(j)-S_lag(:,2))<1e-8);%(S_rov_aut(j)==S_lag(:,2)) ;
                                S_lag_1=(abs(S(j,1)-S_lag(:,1))<1e-8);%(S(j,1)==S_lag(:,1));
                                S_lag_ind=find(S_lag_1.*S_lag_2==1);
                                % for the individual, autarky value is utility of
                                % current consumption plus relevant trans
                                % probability times waut_lag, the vector of values
                                % from the next smaller insurance group
                                if one_per_aut==1
                                    waut(j,1)=utility(rho,caut(j,1)*(1-Pun_fac))+beta*P_lag(S_lag_ind,:)*waut_lag(:,1);
                                    % for rest of village, the autarky value is just
                                    % expected disc utility of current consumption
                                    if Rov_sb==0
                                        waut(j,2)=utility(rho,caut(j,2)*(1-Pun_fac))+beta*P(j,:)*eeewaut(:,2);
                                    else
                                        waut(j,2)=utility(rho,caut(j,2)*(1-Pun_fac))+beta*P(j,:)*waut_rov_sb;
                                    end
                                else
                                    waut(j,1)=waut_lag(S_lag_ind,1);
                                    if Rov_sb==0
                                        waut(j,2)=eeewaut(j,2);
                                    else
                                        waut(j,2)=waut_rov_sb(j,1);
                                    end
                                end
                                
                                for g=1:size(S_rov_aut_all,4)
                                    if STABLE(g)==1
                                        for k=1:size(S_rov_aut_all,3)
                                            S_lag_2=(S_rov_aut_all(j,1,k,g)==S_LAG(:,2,g)) ;
                                            S_lag_1=(S(j,1)==S_LAG(:,1,g));
                                            S_lag_ind=find(S_lag_1.*S_lag_2==1);
                                            
                                            agentstate=[];
                                            for kk=1:length(incomevals)
                                                if idio_dist_aut_all(j,kk,k,g)>0
                                                    agentstate=[agentstate kk];
                                                end
                                            end
                                            
                                            if size(S_lag_ind)>0
                                                for kk=1:length(agentstate) %%%note, you only need one entry per income state, not one entry per agent!!!!!
                                                    
                                                    if one_per_aut==1
                                                        waut_check_ind(j,1,k,g)=utility(rho,S(j,1))+beta*P_LAG(S_lag_ind,find((P_LAG(S_lag_ind,:,g)>0)),g)*WAUT_LAG(find((P_LAG(S_lag_ind,:,g)>0)),1,g);
                                                        waut_check(j,kk,k,g)=utility(rho,incomevals(agentstate(kk)))+beta*P_LAG(S_lag_ind,find((P_LAG(S_lag_ind,:,g)>0)),g)*WAUT_LAG(find((P_LAG(S_lag_ind,:,g)>0)),2,g);
                                                    else
                                                        waut_check_ind(j,1,k,g)=WAUT_LAG(S_lag_ind,1,g);
                                                        waut_check(j,kk,k,g)=WAUT_LAG(S_lag_ind,2,g);
                                                    end
                                                    S_check(j,kk,k,g)=agentstate(kk); 
                                                end
                                            end
                                        end
                                    end
                                end
                            end
                        end
                        
                        % ====================================================
                        % Compute constrained insurance
                        % ====================================================
                        gammagrid=[lambda vss/ScaleRov*(1-lambda)];
                        gammagrid=[gammagrid ones(m,1) gammagrid(:,2)./gammagrid(:,1)];
                        gammagridold=gammagrid;
                        clear gammagrid
                        gammagrid=[linspace(1.2*max(S(:,1)),min(S(:,1)),m)',linspace(min(S(:,2)),1.2*max(S(:,2)),m)',ones(m,1),linspace(1.2*max(S(:,1)),min(S(:,1)),m)'./linspace(min(S(:,2)),1.2*max(S(:,2)),m)'];
                        
                        % ====== make grids of weights, consumption, and Lagrange multipliers for updating ======
                        n=length(gammagrid);
                        w=zeros(n,state,agent);
                        gammar=zeros(n,state,agent);
                        c=zeros(n,state,agent);
                        phi=zeros(n,state,agent);
                        
                        % 'first-best' consumption given aggregate resources and planner weights
                        for j=1:state
                            c(:,j,1)=Y(j)./(1+ScaleRov*gammagrid(:,4));
                            c(:,j,2)=(Y(j)-c(:,j,1))/ScaleRov;
                        end
                        cfirstbest=c;
                        %  relative planner weights
                        for j=1:state
                            for k=1:agent
                                gammar(:,j,k)=gammagrid(:,agent+k);
                            end
                        end
                        % Compute the first best value functions without binding
                        % constraints
                        w=zeros(n,state,agent); dww=1;
                        wold=1/(1-beta)*utility(rho,cfirstbest(:,:,:));
                        while dww>0.00001
                            w(:,:,1)=utility(rho,cfirstbest(:,:,1))+beta*wold(:,:,1)*P';
                            w(:,:,2)=utility(rho,cfirstbest(:,:,2))+beta*wold(:,:,2)*P';
                            dww=max(max(max(abs(w-wold))));
                            wold=w;
                        end
                        wfirstbest=w;
                        % continuation utility as function of planner weight and state of the
                        % economy
                        ewfirstbest(:,:,1)=wfirstbest(:,:,1)*P';
                        ewfirstbest(:,:,2)=wfirstbest(:,:,2)*P';
                        
                        % ====================================================
                        % Start the loop with the enforcement constraints
                        % ====================================================n
                        ewsb=ewfirstbest;
                        wsb=wfirstbest;
                        options=optimset('Display','off');
                        eps=0.0000001;t=0;itgap=0; GAP=2; gap=2; gaphist=gap; gapincrease=0;
                        csolold=cfirstbest;
                        if beta<0.85
                            gapcrit=gapcritset;
                            maxitgap=maxitgapset;
                        else
                            gapcrit=gapcritset*1
                            maxitgap=maxitgapset;
                        end
                        
                        %%idiosyncratic distribution in rov for all states (note this code
                        %%is correct for three individual income states, otherwise have to
                        %%adjust)
                        while GAP >=gapcrit && itgap<=maxitgap
                            rr=[];
                            rr1=[];
                            rr2=[];
                            rrnone=[];
                            rrboth=[];
                            nonconv=zeros(state,3)*nan;
                            itgap=itgap+1;
                            coalnonsust=nan(state,2);
                            count=0; countgrid=[ones(n,state)];
                            for j=1:state
                                
                                %%%for each grid point, figure out which set of
                                %%%individuals from rest of village would want to come
                                %%%along with individual.
                                
                                %%instead of solving problem for each grid point,
                                %%note that there are only three possibilities: 1)
                                %%Constraint binding for Person 1 only, the
                                %%solution is the same for each lambda for which
                                %%this is the case
                                
                                index=[1 j];
                                % 1. constraint binding for agent 1
                                % (individual)
                                rr1=find(utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)<waut(j,1) & utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)>=waut(j,2)) ;
                                if vss>1 & coaldev==1
                                    rr1cd=[];
                                    for g=size(waut_check,4):-1:0 %%%check from the largest group down, 0 is an individual deviation
                                        
                                        if  STABLECHECK(g+1)==1   %%%you need to add one here, because you've got the individual in this vector as well
                                            %%%check the largest possible coalition first:
                                            if g>0
                                                for k=1:size(waut_check,3)   %%%%this code relies on the most binding coalition being with the individual with the highest income, may no longer be true with negative serial correlation
                                                    if waut_check(j,1,k,g)~=0
                                                        coalition=[utility(rho,cfirstbest(rr1,j,1))+beta*ewsb(rr1,j,1)-waut_check_ind(j,1,k,g)];
                                                        for kk=1:size(waut_check,2)
                                                            if waut_check(j,kk,k,g)~=0
                                                                coalition=[coalition utility(rho,cfirstbest(rr1,j,2))+beta*ewsb(rr1,j,2)-waut_check(j,kk,k,g)];
                                                            end
                                                        end
                                                        %%%Step 1: find the grid points for which you are now trying to find a solution
                                                        rrcheck=coalition<0;
                                                        rrcheck=prod(rrcheck,2);
                                                        rrcheck=find(rrcheck>0);
                                                        if size(rrcheck>0)
                                                            bind=1;
                                                           %[j k g waut_check_ind(j,1,k,g)]
                                                            %%%%Step 2: find the solution
                                                            coalitionsolve=[utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)-waut_check_ind(j,1,k,g)];
                                                            for kk=1:size(waut_check,2)
                                                                if waut_check(j,kk,k,g)~=0
                                                                    coalitionsolve=[coalitionsolve utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)-waut_check(j,kk,k,g)];
                                                                end
                                                            end
                                                            rrchecksolve=coalitionsolve<0;
                                                            rrchecksolve=prod(rrchecksolve,2);
                                                            rrchecksolve=find(rrchecksolve>0);
                                                            rr=rrchecksolve(end)+1;
                                                            if rr>length(gammagrid)
                                                                coalnonsust(j,1)=-2;
                                                                output=[nan,nan,-50];
                                                            else
                                                                clear x0
                                                                
                                                                x0(1) = cfirstbest(rr,j,1);  %%redundant because we no longer need an initial guess
                                                                %T vectors of consumption and utilty that meet
                                                                %participation constraints
                                                                [x, output]=agent1maxrho_rov_06_Mar(x0,index,ewsb,[waut_check_ind(:,1,k,g) waut(:,2)],bind,Y,n,beta,rho,ScaleRov,cfirstbest,rr,coaldev);
                                                                
                                                                %This keeps track of nonconvergence, but
                                                                %allows the programme to continue - sometimes it
                                                                %converges in a subsequent iteration
                                                            end
                                                            if output(3) <=0
                                                                nonconv(j,1)=output(3);
                                                                csol(rr1(rrcheck),j,1)=ones(size(rr1(rrcheck),1),1)*caut(j,1);
                                                                csol(rr1(rrcheck),j,2)=ones(size(rr1(rrcheck),1),1)*caut(j,2);
                                                                gammar(rr1(rrcheck),j,2)=(Y(j)-caut(j,1))/caut(j,1)/ScaleRov;
                                                                wsbnewsol(rr1(rrcheck),j,1)=ones(size(rr1(rrcheck),1),1)*waut(j,1);
                                                                wsbnewsol(rr1(rrcheck),j,2)=ones(size(rr1(rrcheck),1),1)*waut(j,2);
                                                                phisol(rr1(rrcheck),j,1)=gammagrid(rr1(rrcheck),4).^rho./gammar(rr1(rrcheck),j,2).^rho-1;
                                                                phisol(rr1(rrcheck),j,2:agent)=0	;
                                                            else
                                                                csol(rr1(rrcheck),j,1)=x(1);
                                                                gammar(rr1(rrcheck),j,2)=(Y(j)-x(1))/x(1)/ScaleRov;
                                                                csol(rr1(rrcheck),j,2)=(Y(j)-x(1))/ScaleRov;
                                                                wsbnewsol(rr1(rrcheck),j,1)=output(1);
                                                                wsbnewsol(rr1(rrcheck),j,2)=output(2);
                                                                phisol(rr1(rrcheck),j,1)=gammagrid(rr1(rrcheck),4).^rho./gammar(rr1(rrcheck),j,2).^rho-1;
                                                                phisol(rr1(rrcheck),j,2:agent)=0	;
                                                            end
                                                        end
                                                    end
                                                end
                                                
                                                rr1(rrcheck)=[];  %%%because these grid points have been dealt with for coalitional deviation of size g
                                                rrnoneadd=rr1; %%%these are the grid points you still need to deal with for coalitional deviations of size smaller g. 
                                            elseif g==0 %%%only an individual deviation
                                                coalition=[utility(rho,cfirstbest(rr1,j,1))+beta*ewsb(rr1,j,1)-eeewaut(j,1)];
                                                %find where there is a binding constraint                                              
                                                rrcheck=coalition<0;
                                                rrcheck=prod(rrcheck,2);
                                                rrcheck=find(rrcheck>0);
                                                if size(rrcheck>0)
                                                    bind=1;
                                                    %find where the
                                                    %solution is
                                                    coalitionsolve=[utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)-eeewaut(j,1)];
                                                    rrchecksolve=coalitionsolve<0;
                                                    rrchecksolve=prod(rrchecksolve,2);
                                                    rrchecksolve=find(rrchecksolve>0);
                                                    
                                                    rr=rrchecksolve(end)+1;
                                                    if rr>length(gammagrid)
                                                        coalnonsust(j,1)=-2;
                                                        output=[nan,nan,-50];
                                                    else
                                                        clear x0
                                                        x0(1) = cfirstbest(rr,j,1);  %%redundant because we no longer need an initial guess
                                                        [x, output]=agent1maxrho_rov_06_Mar(x0,index,ewsb,[eeewaut],bind,Y,n,beta,rho,ScaleRov,cfirstbest,rr,coaldev);
                                                        
                                                        %This keeps track of nonconvergence, but
                                                        %allows the programme to continue - sometimes it
                                                        %converges in a subsequent iteration
                                                    end
                                                    if output(3) <=0
                                                        nonconv(j,1)=output(3);
                                                        csol(rr1(rrcheck),j,1)=ones(size(rr1(rrcheck),1),1)*caut(j,1);
                                                        csol(rr1(rrcheck),j,2)=ones(size(rr1(rrcheck),1),1)*caut(j,2);
                                                        gammar(rr1(rrcheck),j,2)=(Y(j)-caut(j,1))/caut(j,1)/ScaleRov;
                                                        wsbnewsol(rr1(rrcheck),j,1)=ones(size(rr1(rrcheck),1),1)*waut(j,1);
                                                        wsbnewsol(rr1(rrcheck),j,2)=ones(size(rr1(rrcheck),1),1)*waut(j,2);
                                                        phisol(rr1(rrcheck),j,1)=gammagrid(rr1(rrcheck),4).^rho./gammar(rr1(rrcheck),j,2).^rho-1;
                                                        phisol(rr1(rrcheck),j,2:agent)=0	;
                                                    else
                                                        csol(rr1(rrcheck),j,1)=x(1);
                                                        gammar(rr1(rrcheck),j,2)=(Y(j)-x(1))/x(1)/ScaleRov;
                                                        csol(rr1(rrcheck),j,2)=(Y(j)-x(1))/ScaleRov;
                                                        wsbnewsol(rr1(rrcheck),j,1)=output(1);
                                                        wsbnewsol(rr1(rrcheck),j,2)=output(2);
                                                        phisol(rr1(rrcheck),j,1)=gammagrid(rr1(rrcheck),4).^rho./gammar(rr1(rrcheck),j,2).^rho-1;
                                                        phisol(rr1(rrcheck),j,2:agent)=0	;
                                                    end
                                                end
                                                rr1(rrcheck)=[];  %%%because these grid points have been dealt with
                                                rrnoneadd=rr1; % these are additional gridpoints where no one has a binding constraint. 
                                            end
                                        end
                                    end
                                else
                                    %%%%in first group size iteration, i.e. for vss=1
                                    rr1cd=rr1;
                                    rr1id=[];
                                    rrnoneadd=[];
                                end
                                
                                %%%%this is the loop in case CD vss=1 OR
                                %%%%ID
                                if size(rr1cd)>0
                                    bind=1;
                                    rr=rr1cd(end)+1;
                                    if rr>length(gammagrid)
                                        coalnonsust(j,1)=-2;
                                        output=[nan,nan,-50];
                                    else
                                        if size(rr1cd)>0
                                            clear x0
                                            x0(1) = cfirstbest(rr1cd(1),j,1);  %%redundant because we no longer need an initial guess
                                            [x, output]=agent1maxrho_rov_06_Mar(x0,index,ewsb,waut,bind,Y,n,beta,rho,ScaleRov,cfirstbest,rr,coaldev);
                                        end
                                        
                                    end
                                    if output(3) <=0
                                        nonconv(j,1)=output(3);
                                        csol(rr1cd,j,1)=ones(size(rr1cd,1),1)*caut(j,1);
                                        csol(rr1cd,j,2)=ones(size(rr1cd,1),1)*caut(j,2);
                                        gammar(rr1cd,j,2)=(Y(j)-caut(j,1))/caut(j,1)/ScaleRov;
                                        wsbnewsol(rr1cd,j,1)=ones(size(rr1cd,1),1)*waut(j,1);
                                        wsbnewsol(rr1cd,j,2)=ones(size(rr1cd,1),1)*waut(j,2);
                                        phisol(rr1cd,j,1)=gammagrid(rr1cd,4).^rho./gammar(rr1cd,j,2).^rho-1;
                                        phisol(rr1cd,j,2:agent)=0	;
                                    else
                                        csol(rr1cd,j,1)=x(1);
                                        gammar(rr1cd,j,2)=(Y(j)-x(1))/x(1)/ScaleRov;
                                        csol(rr1cd,j,2)=(Y(j)-x(1))/ScaleRov;
                                        wsbnewsol(rr1cd,j,1)=output(1);
                                        wsbnewsol(rr1cd,j,2)=output(2);
                                        phisol(rr1cd,j,1)=gammagrid(rr1cd,4).^rho./gammar(rr1cd,j,2).^rho-1;
                                        phisol(rr1cd,j,2:agent)=0	;
                                    end
                                end
                                %%%%this becomes relevant for vss>1, these are all the left-overs that are not dealt with by joint deviations.
                                if size(rrnoneadd)>0
                                    
                                    csol(rrnoneadd,j,1:agent)=cfirstbest(rrnoneadd,j,1:agent);
                                    gammar(rrnoneadd,j,1:agent)=gammagrid(rrnoneadd,agent+1:2*agent);
                                    phisol(rrnoneadd,j,1:agent)=0;
                                    wsbnewsol(rrnoneadd,j,1)=utility(rho,cfirstbest(rrnoneadd,j,1))+beta*ewsb(rrnoneadd,j,1);
                                    wsbnewsol(rrnoneadd,j,2)=utility(rho,cfirstbest(rrnoneadd,j,2))+beta*ewsb(rrnoneadd,j,2);
                                end
                                
                                % 2. Constraint binding for agent 2 (RoV), I actually
                                % think you don't need to distinguish between whether
                                % it is here binding for agent 1 or not, in either
                                % case, you need to check that the proposed solution
                                % means agent 1 does not have profitable deviations
                                % left, so I am now combining rr2 and rrboth into one
                                % grid point
                                rr2=find(utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)<waut(j,2)) ;
                                
                                if size(rr2)>0
                                    bind=2; %%i.e. agent 2 has the binding constraint
                                    rr=find(utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)>=waut(j,2));
                                    if isempty(rr)==1
                                        coalnonsust(j,2)=-2;
                                        output=[nan,nan,-50];
                                    else
                                        clear x0
                                        x0(1) = cfirstbest(rr(end),j,1); %%redundant because we no longer need an initial guess
                                        [x, output]=agent1maxrho_rov_06_Mar(x0,index,ewsb,[eeewaut(:,1) waut(:,2)],bind,Y,n,beta,rho,ScaleRov,cfirstbest,rr,coaldev);
                                        %%%%check that agent 1 does not have a binding
                                        %%%%constraint at this new solution, and note that
                                        %%%%the rov must never have less than waut(j,2),
                                        %%%%so if individual has a binding constraint that
                                        %%%%means rov would have to decrease its
                                        %%%%consumption, not possible, so you do not need
                                        %%%%to solve, just check!
                                    end
                                    % output(3) equals 1 if constraints
                                    % are satisfied, but -2 if not
                                    if output(3)<=0
                                        nonconv(j,2)=output(3)*100;
                                        csol(rr2,j,1)=ones(size(size(rr2,1),1),1)*caut(j,1);
                                        gammar(rr2,j,2)=(Y(j)-caut(j,1))/caut(j,1)/ScaleRov;
                                        csol(rr2,j,2)=ones(size(size(rr2,1),1),1)*caut(j,2);
                                        wsbnewsol(rr2,j,1)=ones(size(rr2))*waut(j,1);
                                        wsbnewsol(rr2,j,2)=ones(size(rr2))*waut(j,2);
                                        phisol(rr2,j,2)=gammar(rr2,j,2).^rho./gammagrid(rr2,4).^rho-1;
                                        phisol(rr2,j,1)=0;
                                    else
                                        %%%%first check that agent 1 (in conjunction with the rest of the village)
                                        %%%%truly does not have binding constraints
                                        %%%%at the proposed solution
                                        if vss>1  & coaldev==1 % only need to check this once coalitional deviations become relevant
                                            for g=size(waut_check,4):-1:1 %%%check from the largest group down, 0 is an individual deviation
                                                
                                                if  STABLECHECK(g+1)==1   %%%you need to add one here, because you've got the individual in this vector as well
                                                    %%%check the largest possible coalition first:
                                                    for k=1:size(waut_check,3)   %%%%this code relies on the most binding coalition being with the individual with the highest income, may no longer be true with negative serial correlation
                                                        if waut_check(j,1,k,g)~=0
                                                            coalition=[output(1)-waut_check_ind(j,1,k,g)];
                                                            for kk=1:size(waut_check,2)
                                                                if waut_check(j,kk,k,g)~=0
                                                                    coalition=[coalition output(2)-waut_check(j,kk,k,g)];
                                                                end
                                                            end
                                                            %%%Step 1: find the grid points for which you are now trying to find a solution
                                                            rrcheck=coalition<0;
                                                            rrcheck=prod(rrcheck,2);
                                                            rrcheck=find(rrcheck>0);
                                                            if size(rrcheck)>0
                                                                output(3)=-2;
                                                            end
                                                        end
                                                    end
                                                end
                                            end
                                        end
                                        
                                        if output(3)>0
                                            csol(rr2,j,1)=x(1);
                                            gammar(rr2,j,2)=(Y(j)-x(1))/x(1)/ScaleRov;
                                            csol(rr2,j,2)=(Y(j)-x(1))/ScaleRov;
                                            wsbnewsol(rr2,j,1)=output(1);
                                            wsbnewsol(rr2,j,2)=output(2);
                                            phisol(rr2,j,2)=gammar(rr2,j,2).^rho./gammagrid(rr2,4).^rho-1;
                                            phisol(rr2,j,1)=0;
                                        else
                                            nonconv(j,2)=output(3)*100;
                                            csol(rr2,j,1)=ones(size(size(rr2,1),1),1)*caut(j,1);
                                            gammar(rr2,j,2)=(Y(j)-caut(j,1))/caut(j,1)/ScaleRov;
                                            csol(rr2,j,2)=ones(size(size(rr2,1),1),1)*caut(j,2);
                                            wsbnewsol(rr2,j,1)=ones(size(rr2))*waut(j,1);
                                            wsbnewsol(rr2,j,2)=ones(size(rr2))*waut(j,2);
                                            phisol(rr2,j,2)=gammar(rr2,j,2).^rho./gammagrid(rr2,4).^rho-1;
                                            phisol(rr2,j,1)=0;
                                        end
                                    end
                                end
                                
                                % 3. Constraint all slack
                                rrnone=find(utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)>=waut(j,1) & utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)>=waut(j,2)) ;
                                if size(rrnone)>0
                                    csol(rrnone,j,1:agent)=cfirstbest(rrnone,j,1:agent);
                                    gammar(rrnone,j,1:agent)=gammagrid(rrnone,agent+1:2*agent);
                                    phisol(rrnone,j,1:agent)=0;
                                    wsbnewsol(rrnone,j,1)=utility(rho,cfirstbest(rrnone,j,1))+beta*ewsb(rrnone,j,1);
                                    wsbnewsol(rrnone,j,2)=utility(rho,cfirstbest(rrnone,j,2))+beta*ewsb(rrnone,j,2);
                                end
                            end
                            gapold=max(norm(wsbnewsol(:,:,1)-wsb(:,:,1)),norm(wsbnewsol(:,:,2)-wsb(:,:,2)));
                            gap1=max(max(max(abs(invutility(beta,rho,wsbnewsol(:,:,1))-invutility(beta,rho,wsb(:,:,1)))./invutility(beta,rho,wsb(:,:,1)))),max(max(abs(invutility(beta,rho,wsbnewsol(:,:,2))-invutility(beta,rho,wsb(:,:,2)))./invutility(beta,rho,wsb(:,:,2)))))/(1-beta);
                            gap2=max(max(max(abs(csol(:,:,1)-csolold(:,:,1))./csol(:,:,1))),max(max(abs(csol(:,:,2)-csolold(:,:,2))./csol(:,:,2))));
                            gap=max(gap1,gap2);
                            gaphist(itgap+1)=gap;
                            GAP=max(gaphist(max(1,length(gaphist)-convlength):end))
                            disp('gap,gapold')
                            disp([gap,gapold])
                            nonconvind(Qbeta,Qrho,vss,incproc)=mean(nonconv(find(~isnan(nonconv))));
                            Coalnonsust(Qbeta,Qrho,vss,incproc)=mean(coalnonsust(find(~isnan(coalnonsust))));
                        
                            if  (min(min(nonconv))<=0 ||min(min(coalnonsust))<0) & itgap>=3 & coaldev==1
                                itgap
                                disp('Autarky or coalition not sustainable')
                                [nonconv,S]
                                for j=1:state
                                    csol(:,j,1)=caut(j,1);
                                    csol(:,j,2)=caut(j,2);
                                end
                                break
                            end
                            csolold=csol;
                            wsb=wsbnewsol;
                            ewsb(:,:,1)=(P*wsb(:,:,1)')';
                            ewsb(:,:,2)=(P*wsb(:,:,2)')';
                        end
                        CRITGAP(vss)=gaphist(itgap-1);
                        if itgap>=maxitgap
                            nonconvind(Qbeta,Qrho,vss,incproc)=gap;
                        end
                        if gapincrease==1
                            nonconvind(Qbeta,Qrho,vss,incproc)=-gap;
                        end
                        
                        % =============================================================
                        % calculate values of outside option for coalitional deviations
                        % =============================================================
                        
                        if coaldev==1 && (nonconvind(Qbeta,Qrho,vss,incproc)<=0 || Coalnonsust(Qbeta,Qrho,vss,incproc)<=0) && itgap>=3 && vss>1
                            % if vss is unsustainable, the outside option for vss+1 remains that
                            % for vss
                            waut_lag(:,1)=waut_lag(:,1);
                            waut_lag(:,2)=waut_lag(:,2);
                            P_lag=P_lag;
                            S_lag=S_lag;
                            % if the current insurance group is not sustainable, it
                            % cannot be the outside option for the next bigger group,
                            % so increase the 'lagind' indicator by one
                            lagind=lagind+1;
                            gap=-1;
                            STABLE(vss)=0;
                            STABLECHECK=[1 STABLE];
                            % ... or if there is no previous, use just conditional exp of
                            % autarky consumption as calculated above
                        elseif (nonconvind(Qbeta,Qrho,vss,incproc)<=0 || Coalnonsust(Qbeta,Qrho,vss,incproc)<=0)  && itgap>=3 && (vss==1 || coaldev==0)
                            waut_lag(:,1)=eeewaut(:,1);
                            waut_lag(:,2)=eeewaut(:,2);
                            P_lag=P;
                            S_lag=S;
                            lagind=1;
                            gap=-1;
                            
                            S_rov_KEEP(1:state/3,vss)=S_rov;
                            S_KEEP(1:state,:,vss)=S;
                            PHI_KEEP(:,1:state,:,vss)=phisol;
                            STABLE(vss)=1;
                            STABLECHECK=[1 STABLE];
                            IDIO_DIST_ROV_KEEP(1:state/3,1:3,vss)=idio_dist;
                        else
                            lagind=1;
                            % when ins group of size vss is sustainable, its value
                            % becomes outside option for next bigger group size
                            for j=1:state
                                % outside option for next larger insurance group is
                                % value with equal weights vss ...
                                eewaut(j,1)=wsb((n+1)/2,j,1);
                                eewaut(j,2)=wsb((n+1)/2,j,2);
                            end
                            waut_lag=eewaut;
                            S_lag=S;
                            P_lag=P;
                            S_rov_KEEP(1:state/3,vss)=S_rov;
                            S_KEEP(1:state,:,vss)=S;
                            PHI_KEEP(:,1:state,:,vss)=phisol;
                            WAUT_LAG(1:state,1:2,vss)=waut_lag;
                            P_LAG(1:state,1:state,vss)=P_lag;
                            S_LAG(1:state,:,vss)=S_lag;
                            STABLE(vss)=1;
                            STABLECHECK=[1 STABLE]; %%%this includes stability for the individual as well, makes the loop above easier
                            IDIO_DIST_ROV_KEEP(1:state/3,1:3,vss)=idio_dist;
                        end
                        if itgap>=maxitgap
                            nonconvind(Qbeta,Qrho,vss,incproc)=gap;
                        end
                        
                    elseif coaldev==3 % full insurance with coalitional deviations
                        % stationary distribution of aggregate income states
                        SS_S=P^300;
                        SSS_S=SS_S(1,:);
                        SSS_Ssave(vss+1,1:length(SSS_S))=SSS_S;
                        sizeSSS_Ssave(vss+1)=length(SSS_S);
                        % value from consuming average income is V_FI
                        if vss==1
                            sizeSSS_Ssave(1)=3;
                            SSS_Ssave(1,1:3)=SS;
                            V_FI(1,1:3)=(inv(diag(ones(length(incomevals),1))-beta*Q1)*utility(rho,incomevals));
                            P_LAG(1:3,1:3,vss)=Q1; % array with transition matrices for smaller groups
                            Y_lag(1,1:3,1)=incomevals'; % array with income vectors for smaller groups
                        end
                        V_FI(vss+1,1:length(S))=(inv(diag(ones(length(S),1))-beta*P)*utility(rho,Y/(vss+1)))';
                        
                        % vvss are the number of high income individuals
                        Y_lag(vss+1,1:length(Y))=Y;
                        P_LAG(1:length(S),1:length(S),vss+1)=P;
                        % vvss are the number of high income individuals,
                        % required to correspond to stable group sizes. NB:
                        % individuals and pairs are always stable
                        for vvss=[1,stable'+1]
                            ind1=find(abs(Y-(vvss)*max(incomevals)-(vss+1-vvss)*min(incomevals))<0.00000001); % this is the indicator of the state today with vvss high and vss+1-vvss low income individuals - can be several
                            ind2=find(abs(Y_lag(vvss,:)-(vvss)*max(incomevals))<0.000000001); % this is the indicator of the State in a group of vvss individuals associated with only high income individuals
                            % this is the utility loss from making transfers today when
                            % vvss out of vss individuals have high income and vss-vvss
                            % have minimum income
                            deltav(vss,vvss)=(utility(rho,max(incomevals))-utility(rho,(vvss*max(incomevals)+(vss+1-vvss)*min(incomevals))/(vss+1)));
                            % the difference in continuatoin values from having full
                            % insurance with vss vs vvss individuals
                            % NB: use ind1(2), which corresponds to the
                            % individuual having high income
                            check1(vss,vvss)=beta*(P(ind1(end),:)*V_FI(vss+1,1:sizeSSS_Ssave(vss+1))');
                            check2(vss,vvss)=beta*(P(ind1(1),:)*V_FI(vss+1,1:sizeSSS_Ssave(vss+1))');
                            deltaV(vss,vvss)=beta*(P(ind1(end),:)*V_FI(vss+1,1:sizeSSS_Ssave(vss+1))'-P_LAG(ind2(1),1:sizeSSS_Ssave(vvss),vvss)*V_FI(vvss,1:sizeSSS_Ssave(vvss))');
                            % if this is negative for vvss, and the group with vvss is stable, then full insurance with vss individuals
                            % is sustainable against a deviation with vvss individuals
                            Delta(vss,vvss)=deltav(vss,vvss)-deltaV(vss,vvss)
                            
                        end
                        Stable(vss,stable)=(Delta(vss,stable)<0);
                        if sum(Delta(vss,:)<0)==length(stable)+1
                            countstable=countstable+1;
                            stable(countstable,1)=vss;
                        end
                        clear yind  y yy dy aggyy aggy
                    end
                end
            end
            
            %==========================================================
            % Simulate the economy
            % =========================================================
            % This lists all possible combinations of individual income
            % values (IDIO_DIST), and the aggregate states (pair of individual income and ROV-income) associated
            % with every individual income value
            
            if coaldev==1
                stable=(find(STABLE==1));
                vss_set_sim=stable;
            elseif coaldev==0
                vss_set_sim=vss_set;
            elseif coaldev==2
                vss_set_sim=vss_set;
            elseif coaldev==3
                vss_set_sim=stable';
            end
            
            for vss=[vss_set_sim]
                if coaldev<2
                    state=STATE(vss);
                    if rov_average==1
                        ScaleRov=vss;
                    else
                        ScaleRov=1;
                    end
                    S=S_KEEP(1:state,:,vss);
                    S_rov=S_rov_KEEP(1:state/3,:,vss);
                    idio_dist=IDIO_DIST_ROV_KEEP(1:state/3,1:3,vss);
                    phisol=PHI_KEEP(:,1:state,:,vss);
                    IDIO_DIST=[ones(state/3,1)*[1 0 0] + idio_dist; ones(state/3,1)*[0 1 0] + idio_dist; ones(state/3,1)*[0 0 1]+idio_dist];
                    distribution=size(IDIO_DIST,1);
                end
                
                Nvill=ceil((vss_high_inddev+1)/(vss+1))
                var_dyagglog=zeros(Nvill*(vss+1),vss_high);var_dyvillagglog=zeros(Nvill*(vss+1),vss_high);vardaggylog=zeros(Nvill*(vss+1),vss_high);betayc=zeros(2,Nvill*(vss+1));vardcshare=zeros(Nvill*(vss+1),vss_high);vardcshare_dylow=zeros(Nvill*(vss+1),vss_high);vardcshare_dyhigh=zeros(Nvill*(vss+1),vss_high);varclog=zeros(Nvill*(vss+1),vss_high);vardclog=zeros(Nvill*(vss+1),vss_high);varylog=zeros(Nvill*(vss+1),vss_high);vardylog=zeros(Nvill*(vss+1),vss_high);varclog_cond=zeros(Nvill*(vss+1),vss_high);varylog_cond=zeros(Nvill*(vss+1),vss_high);
                varc_yhigh_cond=zeros(Nvill*(vss+1),vss_high);varc_ylow_cond=zeros(Nvill*(vss+1),vss_high);vardclog_cond=zeros(Nvill*(vss+1),vss_high);vardylog_cond=zeros(Nvill*(vss+1),vss_high); varc_yhigh=zeros(Nvill*(vss+1),vss_high);
                vardc_dypos_cond=zeros(Nvill*(vss+1),vss_high);
                vardc_dyneg_cond=zeros(Nvill*(vss+1),vss_high); varc_yhigh=zeros(Nvill*(vss+1),vss_high); varc_ylow=zeros(Nvill*(vss+1),vss_high);vardc_dypos=zeros(Nvill*(vss+1),vss_high); vardc_dyneg=zeros(Nvill*(vss+1),vss_high);
                R2_dcdy=zeros(1,Nvill*(vss+1));betadcdy_agg=zeros(3,Nvill*(vss+1));   betadcdy=zeros(2,vss_high);  betadcdy_high=zeros(1,Nvill*(vss+1));betadcdy_low=zeros(1,Nvill*(vss+1));relcvar=zeros(Nvill*(vss+1),vss_high); reldcvar=zeros(Nvill*(vss+1),vss_high); reldc0var=zeros(Nvill*(vss+1),vss_high); relcvar_cond=zeros(Nvill*(vss+1),vss_high); reldcvar_cond=zeros(Nvill*(vss+1),vss_high);
                beta_cond=zeros(2,Nvill*(vss+1),vss_high); beta_ylow_cond=zeros(4,Nvill*(vss+1),vss_high); beta_yhigh_cond=zeros(4,Nvill*(vss+1),vss_high) ;
                beta_agg=zeros(2,Nvill*(vss+1),vss_high); beta_yhigh_agg=zeros(2,Nvill*(vss+1),vss_high); beta_ylow_agg=zeros(2,Nvill*(vss+1),vss_high);
                
                %=============================================
                % Start simulation
                % ============================================
                % number of insurance groups equals village size
                % vss_high+1 divided by size of ins group vss+1
                
                Npanel=Nvill*(vss+1);
                csim0=nan(Nvill,vss+1,T);
                income0=nan(Nvill,vss+1,T);
                income0ind=nan(Nvill,vss+1,T);
                phisim0=nan(Nvill,vss+1,T);
                Yvill=nan(1,T);
                Yagg=nan(Nvill,T);
                Ssim=nan(Nvill,vss+1,T);
                gammarfun=nan(Nvill,vss+1,T);
                phinew=nan(Nvill,vss+1,T);
                indSjj=nan(T,vss+1);
                disp('simulating')
                disp('Reminder: now simulating for village,coaldev,incproc')
                disp([village coaldev incproc])
                disp('now simulating for parameters')
                disp([Qbeta,Qrho])
                disp('vss')
                disp(vss)
                tic
                %initial state
                t=1;
                % initial values of income and consumption (set =
                % income)
                rng(3,'twister');
                shockinittemp=randi(3,Nvill*(vss+1),1);
                shockinit=reshape(shockinittemp,Nvill,vss+1);
                income0ind(:,:,1)=shockinit;
                income0(:,:,1)=incomevals(income0ind(:,:,1));
                csim0(:,:,1)=income0(:,:,1);
                Yagg(:,t)=sum(income0(:,:,t),2);
                IDIO_DIST_sim=[];
                
                while t<T
                    t=t+1;
                    %%%%simulate next period income taking into account
                    %%%% persistence across time
                    rng(t,'twister');
                    shocktemp=rand(Nvill*(vss+1),1);
                    shock=reshape(shocktemp,Nvill,vss+1);
                    for k=1:vss+1
                        rr=(find(shock(:,k)<=cumQ(income0ind(:,k,t-1),1)));
                        income0ind(rr,k,t)=1;
                        rr=(find(cumQ(income0ind(:,k,t-1),1)<shock(:,k) & cumQ(income0ind(:,k,t-1),2)>=shock(:,k)));
                        income0ind(rr,k,t)=2;
                        rr=(find(cumQ(income0ind(:,k,t-1),2)<shock(:,k) & cumQ(income0ind(:,k,t-1),3)>=shock(:,k)));
                        income0ind(rr,k,t)=3;
                        income0(:,k,t)=incomevals(income0ind(:,k,t))  ;
                    end
                    Yagg(:,t)=sum(income0(:,:,t),2);
                end
                
                if coaldev<2
                    t=1;
                    while t<T
                        t=t+1;
                        gammarfun(:,:,t)=(Yagg(:,t-1)*ones(1,vss+1)-csim0(:,:,t-1))./csim0(:,:,t-1)/ScaleRov;
                        gammarfun(:,:,t)=min(gammarfun(:,:,t),max(gammagrid(:,4))*ones(Nvill,vss+1));
                        for ivill=1:Nvill
                            IDIO_DIST_sim=find(((IDIO_DIST(:,1,1)==(sum(income0ind(ivill,:,t)==1))).*(IDIO_DIST(:,2,1)==(sum(income0ind(ivill,:,t)==2))).*(IDIO_DIST(:,3,1)==(sum(income0ind(ivill,:,t)==3))))==1)';
                            for k=1:vss+1
                                % choose for individual k the aggregate state
                                % that applies
                                indSjj(t,k)=IDIO_DIST_sim(S(IDIO_DIST_sim,1)==income0(ivill,k,t));
                                % check that aggregate and individual income0s are
                                % consistent with the states caculated above
                                if abs(Yagg(ivill,t)-(S(indSjj(t,k),1)+ScaleRov*S(indSjj(t,k),2)))>0.000001
                                    stop
                                end
                            end
                            Ssim(ivill,:,t)=indSjj(t,:);
                        end
                        for k=1:vss+1
                            BLA=interp1(gammagrid(:,4),phisol(:,Ssim(:,k,t),1),gammarfun(:,k,t));
                            phinew(:,k,t)=diag(BLA);
                        end
                        
                        gammaraux=((1+phinew(:,:,t))./(1+phinew(:,1,t)*ones(1,vss+1))).^(1/rho).*csim0(:,:,t-1)./(csim0(:,1,t-1)*ones(1,vss+1));
                        csim0(:,:,t)=(gammaraux./(sum(gammaraux,2)*ones(1,vss+1))).*(Yagg(:,t)*ones(1,vss+1));
                    end
                    AA=[];
                    BB=[];
                    CC=[];
                    for t=1:T
                        AA=[AA reshape(csim0(:,:,t)',Nvill*(vss+1),1)];
                        BB=[BB reshape(income0(:,:,t)',Nvill*(vss+1),1)];
                        CC=[CC reshape(income0ind(:,:,t)',Nvill*(vss+1),1)];
                    end
                    csim0=AA;
                    income0ind=CC;
                    income0=BB;
                elseif coaldev==2
                    
                    AA=[];
                    BB=[];
                    CC=[];
                    for t=1:T
                        BB=[BB reshape(income0(:,:,t)',Nvill*(vss+1),1)];
                        CC=[CC reshape(income0ind(:,:,t)',Nvill*(vss+1),1)];
                    end
                    income0ind=CC;
                    income0=BB;
                    sim_save_prob
                elseif coaldev==3
                    t=0;
                    while t<T
                        t=t+1;
                        csim0(:,:,t)=(Yagg(:,t)*ones(1,vss+1))./(vss+1);
                    end
                    AA=[];
                    BB=[];
                    CC=[];
                    for t=1:T
                        AA=[AA reshape(csim0(:,:,t)',Nvill*(vss+1),1)];
                        BB=[BB reshape(income0(:,:,t)',Nvill*(vss+1),1)];
                        CC=[CC reshape(income0ind(:,:,t)',Nvill*(vss+1),1)];
                    end
                    csim0=AA;
                    income0ind=CC;
                    income0=BB;
                end
                disp('simulating takes')
                tsim(Qbeta,Qrho)=toc;
                disp(tsim(Qbeta,Qrho))
                
                % =========================================================
                % calculate moments of consumption and income growth
                % =========================================================
                csim =[]; cshare_raw =[]; dcshare_raw =[]; cc =[]; ccc =[]; yy =[]; aggyy =[]; aggvillincome=[]; yyy_cond  =[]; dc_log_cond =[]; dy_log_cond =[]; ccc_cond=[];
                rng(3, 'twister'); 
                cerr=randn(size(csim0)); logcsim0_groupmeandev=[]; logcsim_groupmeandev=[]; logincome_groupmeandev=[];
                rng(4, 'twister'); 
                incerr=randn(size(csim0));   Momentsest1=nan(35,N_cerr); ind=[];
                % loop over size of cons measurement error
                for kkk=1:size(csim0,1)
                    tempvar(kkk)=var(log(csim0(kkk,:)));
                end
                Tempvar=mean(tempvar);
                for ic=1:N_cerr+1
                    if N_cerr>0
                        csim=exp(log(csim0)+((ic-1)/N_cerr*maxcerr/(1-(ic-1)/N_cerr*maxcerr)*Tempvar)^0.5*cerr);
                    else
                        csim=exp(log(csim0));
                    end
                    logcsim=log(csim);
                    logcsim0=log(csim0);
                    logcsim0_meandev=logcsim0-ones(size(csim0,1),1)*mean(logcsim0,1);
                    logcsim_meandev=logcsim-ones(size(csim0,1),1)*mean(logcsim,1);
                    
                    % add inc meas error back in
                    income=exp(log(income0)+varnu(incproc)^0.5*incerr);
                    logincome=log(income);
                    logincome0=log(income0);
                    logincome_meandev=logincome-ones(size(csim0,1),1)*mean(logincome,1);
                    income0_meandev=income0-ones(size(csim0,1),1)*mean(income0,1);
                    for ivill=1:Nvill
                        aggvillincome((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=ones(vss+1,1)*sum(income((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
                        logaggvillincome=log(aggvillincome);
                        logcsim0_groupmeandev((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=logcsim0((ivill-1)*(vss+1)+1:ivill*(vss+1),:)-ones(vss+1,1)*mean(logcsim0((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
                        logcsim_groupmeandev((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=logcsim((ivill-1)*(vss+1)+1:ivill*(vss+1),:)-ones(vss+1,1)*mean(logcsim((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
                        logincome_groupmeandev((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=logincome((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)-ones(vss+1,1)*mean(logincome((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
                    end
                    aggincome=sum(income,1);
                    % this chooses log or level consumption for
                    % moments
                    if momentclevinlog==1
                        ccc=logcsim(:,Tinit:T);
                        ccc_cond=logcsim_meandev(:,Tinit:T);
                        ccc_groupcond=logcsim_groupmeandev(:,Tinit:T);
                    else
                        ccc=(csim(:,Tinit:T));
                        ccc_cond=(csim_meandev(:,Tinit:T)); %NB this doesn't actually exist...
                    end
                    yyy=logincome(:,Tinit:T);
                    yyy_cond=logincome_meandev(:,Tinit:T);
                    yyy_groupcond=logincome_groupmeandev(:,Tinit:T);
                    cc=logcsim(:,Tinit:T)-logcsim(:,Tinit-1:T-1);
                    var_dc_t_fe(vss)=var(mean(cc,1));
                    cc_cond=logcsim_meandev(:,Tinit:T)-logcsim_meandev(:,Tinit-1:T-1);
                    cc_groupcond=logcsim_groupmeandev(:,Tinit:T)-logcsim_groupmeandev(:,Tinit-1:T-1);
                    yy=logincome(:,Tinit:T)-logincome(:,Tinit-1:T-1);
                    yy_cond=logincome_meandev(:,Tinit:T)-logincome_meandev(:,Tinit-1:T-1);
                    yy_groupcond=logincome_meandev(:,Tinit:T)-logincome_meandev(:,Tinit-1:T-1);
                    aggyy=log(sum(income(:,Tinit:T),1))-log(sum(income(:,Tinit-1:T-1),1));
                    aggvillyy=log(aggvillincome(:,Tinit:T))-log(aggvillincome(:,Tinit-1:T-1));
                    vardaggylog(1,vss)=var(log(aggincome(1,Tinit:T))-log(aggincome(1,Tinit-1:T-1)));
                    
                    % calculate moments for each individiual
                    for ind=1:size(csim,1)
                        dc_log_cond(ind,:)=(ccc_cond(ind,2:end))-(ccc_cond(ind,1:end-1));
                        dy_log_cond(ind,:)=(yyy_cond(ind,2:end))-(yyy_cond(ind,1:end-1));
                        dc_log_groupcond(ind,:)=(ccc_groupcond(ind,2:end))-(ccc_groupcond(ind,1:end-1));
                        dy_log_groupcond(ind,:)=(yyy_groupcond(ind,2:end))-(yyy_groupcond(ind,1:end-1));
                        
                        % conditioning on income increases, and income
                        % decreases
                        % NB: we allocate no change to decrease!!
                        rrhc =[]; rrlc =[]; rrl =[]; rrh=[];
                        rrhc=find(yy(ind,:)>0);
                        rrlc=find(yy(ind,:)<=0);
                        rrhc_cond=find(yy_cond(ind,:)>0);
                        rrlc_cond=find(yy_cond(ind,:)<=0);
                        rrhc_groupcond=find(yy_groupcond(ind,:)>0);
                        rrlc_groupcond=find(yy_groupcond(ind,:)<=0);
                        vardclog(ind,vss)=var(cc(ind,:));
                        vardylog(ind,vss)=var(yy(ind,:));
                        vardc_dypos(ind,vss)=var(cc(ind,rrhc));
                        vardc_dyneg(ind,vss)=var(cc(ind,rrlc));
                        
                        %%%conditional on village income
                        vardclog_cond(ind,vss)=var(dc_log_cond(ind,:));
                        vardylog_cond(ind,vss)=var(dy_log_cond(ind,:));
                        vardc_dypos_cond(ind,vss)=var(dc_log_cond(ind,dy_log_cond(ind,:)>0));
                        vardc_dyneg_cond(ind,vss)=var(dc_log_cond(ind,dy_log_cond(ind,:)<=0));
                        
                        %%%conditional on group income
                        vardclog_groupcond(ind,vss)=var(dc_log_groupcond(ind,:));
                        vardylog_groupcond(ind,vss)=var(dy_log_groupcond(ind,:));
                        vardc_dypos_groupcond(ind,vss)=var(dc_log_groupcond(ind,dy_log_groupcond(ind,:)>0));
                        vardc_dyneg_groupcond(ind,vss)=var(dc_log_groupcond(ind,dy_log_groupcond(ind,:)<=0));
                        
                        vardc_dypos(ind,vss)=var(cc(ind,rrhc));
                        vardc_dyneg(ind,vss)=var(cc(ind,rrlc));
                        
                        %%% Relative variances
                        reldcvar(ind,vss)=var(cc(ind,rrhc))-var(cc(ind,rrlc));
                        reldcvar_cond(ind,vss)=vardc_dypos_cond(ind,vss)-vardc_dyneg_cond(ind,vss);
                        reldcvar_groupcond(ind,vss)=vardc_dypos_groupcond(ind,vss)-vardc_dyneg_groupcond(ind,vss);
                        
                        %%% THE BETA COEFFICIENTS
                        %1) Unconditional
                        %risk-sharing whole sample
                        covtemp=cov(cc(ind,:),aggyy(1,:));
                        betadcdy_agg(ind,vss)=covtemp(2,1)/var(aggyy(1,:));
                        covtemp=cov(cc(ind,:),yy(ind,:));
                        betadcdy(ind,vss)=covtemp(2,1)/(var(yy(ind,:)));
                        [betadcdytemp,temp1,rtemp,temp2,statstemp]=regress(cc(ind,:)',[ones(1,size(yy,2))',aggyy(1,:)',yy(ind,:)']);
                        rrstat(ind,vss)=rtemp(2);
                        betadcdyalt(ind,vss)=betadcdytemp(2);
                        RR2_dcdy(ind,vss)=statstemp(1);
                        clear rtemp betadcdytemp statstemp
                        [betadcdytemp,temp1,rtemp,temp2,statstemp]=regress(cc_cond(ind,:)',[ones(1,size(yy,2))',yy_cond(ind,:)']);
                        rrstat(ind,vss)=rtemp(2);
                        betadcdyalt1(ind,vss)=betadcdytemp(2);
                        RR2_dcdy_cond(ind,vss)=statstemp(1);
                        clear rtemp betadcdytemp statstemp
                        
                        % for those with income gains versus
                        %those with income losses
                        YY=cc(ind,:)';
                        yy_hc(ind,:)=zeros(1,size(cc,2));
                        yy_hc(ind,rrhc)=yy(ind,rrhc);
                        oneone=ones(1,size(cc,2));
                        ones_hc(ind,:)=zeros(1,size(cc,2));
                        ones_hc(ind,rrhc)=oneone(rrhc);
                        XX=[ones(1,size(cc,2))',ones_hc(ind,:)',yy(ind,:)',yy_hc(ind,:)'];
                        beta_yhigh(:,ind,vss)=regress(YY,XX);
                        betadcdy_low(ind,vss)=beta_yhigh(3,ind,vss);
                        betadcdy_high(ind,vss)=beta_yhigh(4,ind,vss)+beta_yhigh(3,ind,vss);
                        
                        %2) conditional on village income
                        covtemp=[];
                        X=[ones(size(yy_cond(ind,:),2),1),yy_cond(ind,:)'];%
                        YY=cc_cond(ind,:)';
                        beta_cond(:,ind,vss)=inv(X'*X)*(X'*YY);
                        X=[ones(size(aggyy(1,:),2),1),aggyy'];%
                        YY=cc(ind,:)';
                        beta_agg(:,ind,vss)=inv(X'*X)*(X'*YY);
                        %for those with income gains versus
                        %those with income losses
                        YY_cond=cc_cond(ind,:)';
                        yy_hc_cond(ind,:)=zeros(1,size(cc_cond,2));
                        yy_hc_cond(ind,rrhc_cond)=yy_cond(ind,rrhc_cond);
                        oneone=ones(1,size(cc_cond,2));
                        ones_hc_cond(ind,:)=zeros(1,size(cc_cond,2));
                        ones_hc_cond(ind,rrhc_cond)=oneone(rrhc_cond);
                        XX_cond=[ones(1,size(cc_cond,2))',ones_hc_cond(ind,:)',yy_cond(ind,:)',yy_hc_cond(ind,:)'];
                        beta_yhigh_cond(:,ind,vss)=regress(YY_cond,XX_cond);
                        
                        %3) conditional on group income
                        covtemp=[];
                        X=[ones(size(yy_groupcond(ind,:),2),1),yy_groupcond(ind,:)'];%
                        YY=cc_groupcond(ind,:)';
                        beta_groupcond(:,ind,vss)=inv(X'*X)*(X'*YY);
                        X=[ones(size(aggvillyy(ind,:),2),1),aggvillyy(ind,:)'];%
                        YY=cc(ind,:)';
                        beta_groupagg(:,ind,vss)=inv(X'*X)*(X'*YY);
                        %risk-sharing for those with income gains versus
                        %those with income losses
                        YY_groupcond=cc_groupcond(ind,:)';
                        yy_hc_groupcond(ind,:)=zeros(1,size(cc_groupcond,2));
                        yy_hc_groupcond(ind,rrhc_groupcond)=yy_groupcond(ind,rrhc_groupcond);
                        oneone=ones(1,size(cc_groupcond,2));
                        ones_hc_groupcond(ind,:)=zeros(1,size(cc_groupcond,2));
                        ones_hc_groupcond(ind,rrhc_groupcond)=oneone(rrhc_groupcond);
                        XX_groupcond=[ones(1,size(cc_groupcond,2))',ones_hc_groupcond(ind,:)',yy_groupcond(ind,:)',yy_hc_groupcond(ind,:)'];
                        beta_yhigh_groupcond(:,ind,vss)=regress(YY_groupcond,XX_groupcond);
                    end
                    var_dyagglog(1,vss)=var(aggyy);
                    R2_dcdy=mean(RR2_dcdy_cond(1:size(csim,1),vss));
                    
                    % =====================
                    % Matrix of moments
                    % =====================
                    % col 1-5 are parameters
                    Momentsest1(1:5,ic)=[vss,beta,rho,rho1vec(incproc),varnu(incproc)]';
                    % 6-11 relative variance dc/dy, dc(dy>0)/dc(dy<0);  unconditional (6:7), conditional on village income
                    % (8:9), conditonal on group income (10:11)
                    Momentsest1(6:11,ic)=[mean(vardclog(1:size(csim,1),vss))/mean(vardylog(1:size(csim,1),vss)), mean(reldcvar(1:size(csim,1),vss))...
                        mean(vardclog_cond(1:size(csim,1),vss))/mean(vardylog_cond(1:size(csim,1),vss)), mean(reldcvar_cond(1:size(csim,1),vss))...
                        mean(vardclog_groupcond(1:size(csim,1),vss))/mean(vardylog_groupcond(1:size(csim,1),vss)), mean(reldcvar_groupcond(1:size(csim,1),vss))]';
                    % 12-14: beta, beta high, beta low unconditional
                    Momentsest1(12:14,ic)=[mean(betadcdy(1:size(csim,1),vss)),mean(betadcdy_high(1:size(csim,1),vss)),mean(betadcdy_low(1:size(csim,1),vss))]';
                    % 13_16: beta, beta high, beta low conditional on
                    % village income
                    Momentsest1(15:17,ic)=[mean(beta_cond(2,1:size(csim,1),vss),2),mean(beta_yhigh_cond(3,1:size(csim,1),vss),2)+mean(beta_yhigh_cond(4,1:size(csim,1),vss),2),mean(beta_yhigh_cond(3,1:size(csim,1),vss),2)]';
                    % 18:20: beta, beta high, beta low conditional on
                    % group income
                    Momentsest1(18:20,ic)=[mean(beta_groupcond(2,1:size(csim,1),vss),2),mean(beta_yhigh_groupcond(3,1:size(csim,1),vss),2)+mean(beta_yhigh_groupcond(4,1:size(csim,1),vss),2),mean(beta_yhigh_groupcond(3,1:size(csim,1),vss),2)]';
                    %21:22 are variances of dclog, dylog
                    Momentsest1(21:22,ic)=[mean(vardclog(1:size(csim,1),vss)),mean(vardylog(1:size(csim,1),vss))];
                    %23:28, relative variance dc/dy (dy>0) and var dc/dy(dy<0);
                    %unconditinoal 23:24, conditional on village income 25:26,
                    %conditional on group income 27:28
                    Momentsest1(23:24,ic)=[mean(vardc_dypos(1:size(csim,1),vss));mean( vardc_dyneg(1:size(csim,1),vss))];
                    Momentsest1(25:26,ic)=[mean(vardc_dypos_cond(1:size(csim,1),vss));mean(  vardc_dyneg_cond(1:size(csim,1),vss))];
                    Momentsest1(27:28,ic)=[mean(vardc_dypos_groupcond(1:size(csim,1),vss));mean( vardc_dyneg_groupcond(1:size(csim,1),vss))];
                    Momentsest1(30,ic)=[mean( vardylog_cond(1:size(csim,1),vss))];
                    Momentsest1(31,ic)=mean(RR2_dcdy(1:size(csim,1),vss));
                    Momentsest1(32,ic)=mean(vardclog(1:size(csim,1),vss));
                    Momentsest1(33,ic)=mean(vardclog_cond(1:size(csim,1),vss));
                    Momentsest1(34,ic)=var_dc_t_fe(vss);
                    Momentsest1(35,ic)=mean(RR2_dcdy_cond(1:size(csim,1),vss));
                    %%% and some other model-dependent measures
                    if coaldev<2
                        Momentsest1(29,ic)=CRITGAP(vss);
                    elseif coaldev==2
                        Momentsest1(29,ic)=phi;
                    elseif coaldev==3
                        Momentsest1(29,ic)=nan;
                    end
                end
                Momentsest(:,:,Qbeta,Qrho,vss)=[Momentsest1];
                toc
            end
            cd(outputpath)
            save(filetosaveest,'Momentsest');
            cd(programpath)
        end
    end
    
    cd(outputpath)
    if rhogrid==0 && buildup==0 && Rov_sb==0 && scatter == 0
        eval(['save ', filetosave, int2str(village),int2str(coaldev), int2str(incproc),int2str(maxincerr==0),int2str(unconditional) int2str(fixedeffect) int2str(one_per_aut) '.mat'])
    elseif Rov_sb==1
        eval(['save SB', filetosave, int2str(village),int2str(coaldev), int2str(incproc),int2str(maxincerr==0),int2str(unconditional) int2str(fixedeffect) int2str(one_per_aut) '.mat'])
    elseif rhogrid==1
        eval(['save Grid', filetosave, int2str(village),int2str(coaldev), int2str(incproc),int2str(maxincerr==0),int2str(unconditional) int2str(fixedeffect) int2str(one_per_aut) '.mat'])
    elseif buildup==1
        eval(['save Buildup', filetosave, int2str(village),int2str(coaldev), int2str(incproc),int2str(maxincerr==0), int2str(unconditional) int2str(fixedeffect) int2str(one_per_aut) '.mat'])
    elseif scatter==1
        eval(['save Scatter', filetosave, int2str(village),int2str(coaldev), int2str(incproc),int2str(maxincerr==0), int2str(unconditional) int2str(fixedeffect) int2str(one_per_aut) '.mat'])
    end
    cd(programpath)
end